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ABSTRACT 

We present new optical long-slit data along 6 position angles of the bulge region of M31. We derive accurate stellar and gas kinematics 
reaching 5 arcmin from the center, where the disk light contribution is always less than 30%, and out to 8 arcmin along the major axis, 
where the disk makes 55% of the total light. We show that the velocity dispersions of McElroy (1983) are severely underestimated 
(by up to 50 km/s). As a consequence, previous dynamical models have underestimated the stellar mass of M3rs bulge by a factor 
2. As a further consequence, the light-weighted velocity dispersion of the galaxy grows to 166 km/s and to 170 km/s if also rotation 
is taken into account, thus reducing the discrepancy between the predicted and measured mass of the black hole at the center of M31 
from a factor 3 to a factor 2. The kinematic position angle varies with distance, pointing to triaxiality, but a quantitative conclusion 
can be reached only after simultaneous proper dynamical modeling of the bulge and disk components is performed. We detect gas 
counterrotation near the bulge minor axis. We measure eight emission-corrected Lick indices. They are approximately constant on 
circles. Using simple stellar population models we derive the age, metallicity and cf-element overabundance profiles. Except for the 
region in the inner arcsecs of the galaxy, the bulge of M31 is uniformly old (> 12 Gyr, with many best-fit ages at the model grid 
limit of 15 Gyr), slightly a-elements overabundant ([a/Fe] a 0.2) and at solar metallicity, in agreement with studies of the resolved 
stellar components. The predicted u-g, g-r and r-i Sloan color profiles match reasonably well the dust-corrected observations, within 
the known limitations of current simple stellar population models. The stellar populations have approximately radially constant mass- 
to-light ratios (M/Lr « 4 - A.SMq/Lq for a Kroupa IMF), in agreement with stellar dynamical estimates based on our new velocity 
dispersions. In the inner arcsecs the luminosity-weighted age drops to 4-8 Gyr, while the metallicity increases to above 3 times 
the solar value. Starting from 6 arcmin from the center along the major axis, the mean age drops to < 8 Gyr, with slight supersolar 
metallicity (« +0. 1 dex) and a-element overabundance (« -1-0.2 dex), for a mass-to-light ratio M/Ljt < 3Mq/Lq. Diagnostic diagrams 
based on the [0III]/Hy6 and [NI]/Hy6 emission line equivalent widths (EWs) ratios indicate that the gas is ionized by shocks outside 10 
arcsec, but an AGN-like ionizing source could be present near the center. We speculate that a gas-rich minor merger happened some 
100 Myr ago, causing the observed minor axis gas counterrotation, the recent star formation event, and possibly some nuclear activity. 

Key words. Galaxies: individual: M31, kinematics and dynamics, stellar content 

1. Introduction try and kinematics is still missing. A bar could also be present 

(Athanassoula & Beaton |2006' Beaton et al. '20071). Moreover, 

This is the first of two papers presenting new optical spectra for investigations of the stellar populations of the central regions of 

the bulge of M31 to study its stellar populations and assess its M31 through the measurement of Lick indices have been per- 

triaxiality through dynamical modeling. Here we present the new formed (Davidge 1997 ). They indicate the presence of a young 

data and constrain the stellar populations. and metal rich population in the inner arcsecs of the galaxy. 

In the past 50 years papers studying the dynamics of our Parallely, studies of the resolved stellar population of the bulge 

neighbour galaxy M31 have been published on a regular ba- of M31 have assessed that the global stellar population of the 

sis, discussing gas kinematics, both by optical spectroscopy M3 1 bulge must be as old as the bulge of Milky Way, resolving 

(Boul esteix et al. ,1987. Pe llet ,1976) . and in HI (Kent 1989a. previous claim s of yo unger ages as due to crowding problems 

Braun [19911 Chemin et al. |2009l and references therein), stel- (Stephens et al. 120031 ). 

lar kinematics concentrating on the central regions to probe the Two considerations convinced us of the necessity to collect 

black hole dynamics (Bender et al. 2005) or considering the new optical spectroscopic information for the bulge of M31, in 

whole bulge (McEkoy 1983 ). The data are used to construct dy- addition to the old age of the dataset of McElroy ( 1983 ). The 

namical models of the galaxy (Widrow et al. [20031 Klypin. Zhao first one is the start of PAndromeda, an extensive monitoring 

and SomerviIIe l2002l ) and possibly probe the tridimensional dis- campaign of M3 1 with the PanSTARRS- 1 telescope and camera 

tribution of its stellar components. The question of the triaxial- system (Kaiser [20041 ). that is in principle able to deliver hun- 

ity of M31 bulge has been posed early on (Stark [T9771 Gerhard dreds of pixel lensing events, probing both bulge and disk re- 

119861 ) and is of significant importance for the understanding of gions. Detailed stellar population and dynamical models, based 

M31, but a definitive quantitative modeling of both photome- on accurate spectral information, are needed to interpret these 
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events as due to a compact baryonic dark matter component (the 
so-called MACHOs) rather than self-lensing of stellar popula- 
tions (Kerins et al. 2001 RifFeser et al. 2006 ). The second is the 
development of new modeling techniques of both simple stellar 
populations and stellar dynamical systems. On the one hand, the 
interpretation of Lick indices (Worthey et al |1994| l in terms of the 
most recent simple stellar population models (Maraston 119981 
I2005I I that take into account the variation of [a/Fe] (Thomas, 
Maraston and Bender 2003 ), allows the accurate determination 
of the stellar population ages, metallicities and overabundances, 
and therefore the prediction of stellar mass-to-light ratios. On the 
other hand, new dynamical modeling codes, like N-MAGIC (De 
Lorenzi et al. 2007 ) allow the flexible dynamical modeling of tri- 
axial structures, optimally exploiting the information contained 
in the line-of-sight velocity distributions that modern programs 
for the analysis of the galaxy optical spectra are able to extract 
(Bender, Saglia and Gerhard 1 1994| l. well beyond the mean ve- 
locities and velocity dispersions of McEkoy (1983J . 

In the following we discuss our new spectroscopic observa- 
tions of the bulge of M31. A future paper (Morganti et al. |2009l l 
will report on the dynamical modeling. In Sect. |2] we present the 
observations and the data reduction. In Sec.[3]we derive the stel- 
lar and gas kinematics and the strengths of the absorption and 
emission lines. In Sec. |4] we discuss the modeling of the new 
spectroscopic data. We analyze the stellar population in Sect. 
14. H and discuss previous axisymmetric dynamical models of the 
bulge of M3 1 in Sect. 14.21 Sect. 14.31 considers the possible ex- 
citation sources compatible with the observed emission line EW 
ratios. We draw our conclusions in Sect. |5] 

2. Observations and Data reduction 

We observed the bulge of M31 using the Low Resolution 
Spectrograph (LRS, Hill et al. [T998] ) at the Hobby-Eberly 
Telescope equipped with a 1.5 arcsec wide, 3.5 arcmin long 
slit, the E2 grism and a Ford Aerospace CCD device, with 
3072x1024 ISfim pixels (usable range 2750x900 pixels) and a 
scale of 0.235 arcsec per pixel. We covered the wavelength range 
/I = 4787 -5857 A with 0.36 A per pixel and an instrumental res- 
olution of cr,„„ - 57 km/s. During the period 10th-19th August 
2007 (see Table [T]i M3 1 was observed in service mode along 6 
position angles. The seeing varies from 1.3 to 2.5 arcsec. At each 
position angle three 10 minutes exposures were taken, the first 
one centered on the galaxy, the second one shifted to the west 
(decreasing RA) 3.5 arcmin in the slit direction, and the third 
one shifted to the east (increasing RA) 3.5 arcmin in the slit di- 
rection. Finally, on the 15th of September 2009 we collected a 
20 minutes exposure along the major axis, shifted 5 arcmin from 
the center in the eastern direction. Fig. [1] shows the distribution 
of the 19 sUt M31 pointings on the sky with the naming conven- 
tion given in Table [T] where MJ shortens for the major axis and 
MN for the minor axis. The position angle of the bulge MJ axis 
PA - 48 is fixed from our 2MASS photometi-y (see Sect.|4!2]and 
Fig- Ell. Note that this is different from the position angle of the 
major axis of the disk of M31 (PA = 38, de Vaucoulers 1 19581 1. 
Apart from the MN pointings, the slits with a given position an- 
gle are always well aligned and with small overlaps, providing 
spectra out to 5 arcmin from the center. In addition, the MJEE 
slit probes the major axis out to 8 arcmin from the center. The 
MNE and MNW are slighly shifted from the MNC central slit 
(by 34" orthogonal to the slit to the west and to the east respec- 
tively). Furthermore, 10 minutes exposures of empty sky regions 
were also taken, as well as several kinematic and Lick standard 
stars, wiggled and trailed along the slit. 



Table 1. Log of the observations. 



Date 


Name 


PA 


FWHM 


Comment 






(deg.) 


(") 




20070812 


MJC 


48 


1.40 


Centered major axis 


20070812 


MJE 


48 


1.46 


Eastern major axis 


20070812 


MJW 


48 


1.38 


Western major axis 


20090915 


MJEE 


48 


2.13 


Eastern major axis 


20090915 


SKYMJEE 


48 


- 


Empty sky 


20070813 


P30C 


78 


1.59 


Centered -1-30° 


20070813 


P30E 


78 


1.91 


Eastern +30° 


20070813 


P30W 


78 


2.08 


Western +30° 


20070813 


P30SKY 


90 


- 


Empty sky 


20070814 


MNC 


138 


1.28 


Centered minor axis 


20070814 


MNE 


138 


1.36 


Eastern minor axis, shifted 


20070814 


MNW 


138 


1.41 


Western minor axis, shifted 


20070814 


MNSKY 


90 


- 


Empty sky 


20070816 


P60C 


108 


1.61 


Centered +60° 


20070816 


P60E 


108 


1.57 


Eastern +60° 


20070816 


P60W 


108 


1.63 


Western +60° 


20070816 


P60SKY 


90 




Empty sky 


20070817 


M30C 


18 


2.53 


Centered -30° 


20070817 


M30E 


18 


2.21 


Eastern -30° 


20070817 


M30W 


18 


2.21 


Western -30° 


20070819 


M60C 


168 


1.43 


Centered -60° 


20070819 


M60E 


168 


1.78 


Eastern -60° 


20070819 


M60W 


168 


1.94 


Western -60° 


20070819 


M60SKY 


90 




Empty sky 




Fig. 1. The positions of the 18 M31 slit pointings superimposed 
on a NOAO image of the Andromeda galaxy (credit Adma 
Block/NOAO/AURA/NSF). See Table[I]for the naming conven- 



The data reduction followed the usual procedure described in 
Mehlert et al. (2000) under MIDAS with some additional steps 
needed to correct for the spectral alignment and the anamor- 
phism of the LSR spectrograph. After bias subtraction and flat- 
fielding, the 2° tilt between the wavelength direction and the 
CCD was removed through appropriate rebinning. Two bad 
columns at positions corresponding to A ^ 4850A were cor- 
rected by interpolation. The wavelength calibration used Cd and 
Ne calibration lamps frames taken at the end of each night, and a 
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3rd order polynomial fitting 9 lines and giving an rms of less than 
one pixel. Residual overall shifts in the starting wavelength of 
the wavelength calibrated science frames were corrected by ref- 
erencing to the 5577 A sky line. After filtering cosmic ray hits, 
the science frames were rebinned to a logarithmic wavelength 
scale. The sky subtraction for the M31 frames was performed 
through the four empty sky frames available. The averaged sky 
continuum amounted to 3.5 counts per pixels (varying by at most 
10% from frame to frame). This is only 14% of the averaged flux 
measured at the extreme ends of the MJE and MJW slits, 23% 
at the end of the MJEE slit, but 44% at those of the MNE and 
MNW slits. In contrast, the flux per pixel in the strongest sky line 
at 5577 A is 160 counts with variations to up to 50% from frame 
to frame. We constructed a line-free and a line-only average sky 
spectrum for each of the available sky frames and subtracted 
it with proper scaling from the corresponding science frames. 
The average of the available sky frames was used for the nights 
without empty sky observations. Stars were extracted tracing the 
peaks corresponding to their trails across the slit as a function 
of wavelength and averaging over a 10 pixel wide window. The 
final spectra are averages of the extracted ones. This procedure 
also allowed us to map the anamorphic distortions of the spec- 
trograph, that at the short and long wavelengths bend the spectra 
upwards in the upper end and downwards in the lower end of 
the slit (see Fig.|2]i. Moreover, the sky was estimated at the ends 
of the slit and subtracted. Finally, the anamorphic mapping was 
used to rebin the sky subtracted frames of M3 1 on a cartesian 
grid. Finally, the M3 1 spectra were rebinned in the radial direc- 
tion to give approximately constant signal-to-noise ratio. 



3. Kinematics and Line Strengtlis 

The stellar kinematics were derived using the Fourier 
Correlation Quotient technique (Bender I1990D , with the addi- 
tion of an iterative procedure to correct for the emission clearly 
detected in the H/3 line and the [OIII] and [NI] doublets, similar 
to the method of Sarzi et al. (I2006I I. The best fitting stellar tem- 
plate convolved to the derived line-of-sight velocity distribution 
(LOSVD) is subtracted from the galaxy spectra. Five Gaussian 
functions are fitted to the residuum spectrum, in correspondence 
to the two doublets and the Hfi line. Six parameters are de- 
rived: the recessional velocity of the gas emission and its veloc- 
ity dispersion (determined by fitting just the [OIII] lines), plus 5 
normalization factors. The best-fitting Gaussians are subtracted 
from the original spectrum and the stellar kinematics is derived 
a second time, using the wavelength range 4870 - 5470 A. The 
final parameters of the gas emission are derived by fitting the 
original spectrum where the improved, iterated best-fitting stel- 
lar kinematic template is subtracted. The best-fitting template 
is chosen from a pool of single star spectra and simple stellar 
population models from Vazdekis ( 119991 ), by minimizing the 
between the galaxy and the LOS VD-convolved template, and re- 
sults in an almost mismatch-free description of the galaxy spec- 
trum. An 8th order polynomial is used to subtract the contin- 
uum and the first and last three channels in Fourier space are 
filtered out. Monte Carlo simulations as described in Mehlert 
et al. ( I2OOOI 1 are performed to derive the final errors on the es- 
timated parameters and estimate the size of systematic errors. 
With the choice of continuum fitting and filtering given above 
the systematic errors are always smaller than the statistical ones. 
We have tested that none of the results presented below depend 
on the precise order of the continuum subtraction within the er- 
rors, as soon as it is larger than 6. The measured kinematics are 



very precise, with statistical errors on velocities of the order of 2 
km/s, of the order of 3 km/s (i.e. 3%), on cr, less than 0.02 on Hj, 
and H4. We investigated the systematic effects of sky subtraction 
by deriving the kinematics from frames where ±10% of the sky 
was additionally subtracted. At the extreme ends of the slits dif- 
ferences in velocities up to 5 km/s, in cr up to 12 km/s, in i/3 up 
to 0.03 and in H4 up to 0.05 are seen. 

Fig. [3] shows the derived stellar line of sight recessional 
velocity V, velocity dispersion cr and Hermite-Gauss coeffi- 
cients Hi and H4 following Bender, Saglia and Gerhard (|1994| l. 
Overall, the kinematic data are symmetric with respect to the 
center within the errors. Stellar rotation is detected at some level 
at all position angles. Along the major axis (MJ) it keeps rising 
with increasing distance from the center, reaching « 100 km/s at 
the last measured point at a; 8 arcmin, where the disk light dom- 
inates. The velocity dispersion reaches a maximum of 170 km/s 
at a; 60 arcsec from the center Along or near the major axis, 
where stellar rotation is clearly detected, the parameter anti- 
correlates (as expected, see Bender et al. |I994| l with the stellar 
velocity out to a 100 arcsec, to change sign at larger distances 
(see discussion in Sect. 14.21 ). In the inner arcsecs the signature 
of the central supermassive black hole is clearly visible, with a 
steep increase in velocity and velocity dispersion. In the inner 
region, the comparison with the kinematics of van der Marel et 
al. (1994) along the major and minor axis shows overall good 
agreement; however they measure systematically (a 0.05) lower 
values of H4. Since we agree with the data of Kormendy and 
Bender ( 1999) and we do not expect residual systematics from 
our Monte Carlo tests (see above), we trust our H4 profiles as 
the correct ones. The comparison with McElroy (1983) shows 
that this old dataset gives seriously underestimated velocity dis- 
persions. The smallest differences are observed along the mi- 
nor axis. The discrepancy is probably due to the sensitivity of 
McElroy analysis's method to the disk light contamination (see 
discussion in Sect. 14. 2b . The recessional velocities are in reason- 
able agreement. 

As a consequence, the mean velocity dispersion used to pre- 
dict the mass of the black hole at the center of M31 (160 km/s, 
Gultekin et al. 120091 corresponding to Mbh - 0.5 x IO^Mq, 
3 times smaller than the measured value) has to be revised 
upwards. Fig. [H top, shows that the (angle averaged) light- 
weighted velocity dispersion converges to 166 km/s and to al- 
most 170 km/s if the rotational velocity is added in quadrature. 
This results in a black hole mass only a factor 2 smaller than ob- 
served. The same figure, bottom, shows that lower values are de- 
rived as a function of the distance from the center, if the weight- 
ing is performed using the total isophotal light. 

The gas kinematics (V^as and cr^ns) are shown in Fig.|5] Gas 
velocities are almost a factor two higher than stellar velocities, 
but less axisymmetric and regular Gas velocity dispersions are 
slightly larger than the instrumental resolution, indicating intrin- 
sic velocity dispersions less than 80 km/s. In the inner 10-20 
arcsec near the minor axis gas counterrotation (i.e., the orienta- 
tion of the gas rotation changes) is observed. The gas velocities 
compare reasonably well with the Ha measurements of Pellet 
(T9W) or Boulesteix et al. ( fT987t . 

We measured Lick line strength index profiles on the M3 1 
emission-corrected spectra as in Mehlert et al. ( 120001 ). from H/3 
to Fe5789. 

As a first check, we verify that we measure line indices uni- 
formly well across the slit. For this purpose, we analyse the spec- 
tra of the Lick standard star HD 72324 (= HR 3369) that was 
observed trailed along the slit. Fig. |6] shows the variation of its 
line indices as a function of the position on the slit. The mea- 
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Fig. 2. The anamorphic distortion of the LRS as mapped by tracing the peaks of the stellar spectra. The traces of 9 stellar spectra 
are shown as full black lines and are compared to the corresponding straight lines (dashed). The interpolated differences (in pixel) 
between the stellar traces and the Cartesian grid as a function of the position on the chip are shown coded by colour 



sured molecular indices Mgl and Mg2 with their widely spread 
continuum windows are increasingly biased towards the ends 
of the slit, probably because of inaccurate spectral flat fielding 
where vignetting becomes important. In contrast, all the other 
indices are well determined and do not vary much . In the fol- 
lowing we do not consider Mgl and Mg2 anymore. As a conse- 
quence, all (atomic) indices presented here below are measured 
in Angstrom. 

As a second step, we assess how well we are on the Lick 
system, since our spectra are not flux calibrated. Our indices 
are measured following the bands definitions of Worthey et al. 
( 1199 4) with the corrections of Trager et al. ( 1998). Fig.Qdemon- 
strates that no tilt nor shift are needed within the errors, with a 
typical rms of 0.2 A. The largest shift, 0.32 A, significant at the 
2.5 sigma level, is observed for the H/3 line, where, as discussed 
above, we have interpolated two bad columns. We have added 
in quadrature the rms of each relation to the errors of the galaxy 
indices. 

Finally, Fig. [8] shows the typical effect of the correction for 
gas emission applied to the spectra. The plot shows the differ- 
ences Al = Icor - Iraw between the emission corrected index 
Icor and the values Iraw measured before the correction. Before 
the correction, the H/? and Fe5015 indices are smaller (since the 
emission is present in the index window) by typically 0.29 A 
and 0.41 A respectively, and Mgb is larger (since the emission 
is present in the redder continuum window) by 0.07 A. While 
the AH/3 and A[NI] match fully the EW of the emission, A[OIII] 
amounts to the difference between the EWs of the [OIII] doublet 
redder line (that falls in the band window of the Fe5015 index) 
and the bluer line (that falls in the bluer continuum window of 
the index). The other indices are not affected. Again, we add in 
quadrature the errors on the emission EWs of Hy6, [OIII] and 
[NI] to the ones of the absorption line strengths H/3, Fe5015 and 
Mgb, respectively. 

Fig. |9] shows the measured Lick indices profiles for M31. 
Apart from the central regions, where rapid variations are mea- 



sured, the indices show very mild gradients, indicative of homo- 
geneous stellar populations (see discussion in Sect. 14. It . Along 
the major axis and at distances larger than k 5 arcmin, where 
the disk light starts to dominate, making 55% of the total at 500 
arcsec from the center, the strenght of the Hp line increases, in- 
dicative of a younger stellar population (see discussion in Sect. 
14. It . Overall, the indices appear not only symmetric with respect 
to the center, but also approximately symmetric on circles. As 
an example, if we compute the rms of the differences of the Mgb 
index measured at the same distance on the major axis and on 
one of the other slit positions, we get values between 0.15 and 
0. 18 A. If we now repeat the procedure interpolating at the same 
isophotal distance (using the ellipticity profile derived in Sect. 
14.21 see Fig. [TsT l. the rms increases to 0.20-0.25 A. The same 
applies to the HjS or the Iron indices. 

Davidge (|1997| l measured Lick indices (not corrected for gas 
emission) in a 22x60 arcsec area centered on M3L His values 
along the centered East- West slit are shown in Fig. [8] on the plots 
for the P30 slit (nearest to the EW position). The agreement is 
good. 

In Fig. [To] we plot the equivalent widths (EWs) of the Hy6 
emission line in A, together with the EW ratios [OIII]/Hy6 and 
[NI]/H/3 of the emission lines considered above in logarithmic 
units. Note that we expect the EW ratios to vary by no more than 
10% with respect to the ratios of the emission line fluxes, due 
to the restricted wavelength range spanned by the three lines. 
Given the logarithmic units, this does not affect the discussion 
given in Sect. 14.31 While the Hy6 and [OIII] lines are almost al- 
ways well detected, with a signal to noise ratio larger than 4, 
the [NI] doublets are often weak, if not undetected. The H/3 EW 
is approximately constant with radius. At radii larger than k 10 
arcsec the same is true for the [OIII]/H/3 ratio (at the value of 3) 
and the [NI]/HyS ratio (at low values smaller than » 0.1). Near the 
center [OIII]/H/3 increases to » 10, combined with an increased 
[NI]/H/3 ratio (^ 1). 
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Fig. 3. The stellar kinematics along different position angles, folded with respect to the center (at velocity -333.2 km/s), antisym- 
metrically for recessional velocities and coefficients, symmetrically for velocity dispersions and H4 coefficients. Filled circles 
refer to radii with increasing RA (East), stars to radii with decreasing RA (West). The number in the velocity panel indicates the 
correction in km/s to the central velocity (-333.2 km/s) applied to achieve maximal symmetry. The plus or minus sign indicates 
whether the velocities of the sides with increasing RA (the filled dots) are receding or approaching. The vertical dotted lines mark 
the transition from the central slit data to those from the outer two. For the major axis, a further vertical dotted line marks the outer 
MJEE dataset. The short-dashed vertical lines mark the seeing values of the central slits (see Table [T]l. The data of McElroy (119831 
at PA=75, 45, 15, 165, 135 and 105 respectively) are indicated by green open circles and triangles. The data of van der Marel et al. 
( 119941 MJ at PA=55, MN at PA=128) are shown by the red continuous lines. The data of Kormendy and Bender (1999) are shown 
in yellow (from their Calcium Triplet spectra) and magenta (from their HyS-Mg-Fe spectra). The blue dashed lines show Model Kl 
of Widrow et al. (i2003) , the blue full lines their Model A (see discussion in Sect. |4.2| l. 



Tables |2]|3] and|4]give format examples of the measured stel- 
lar and gaseous kinematics, and Lick indices as a function of dis- 
tance and position angle, respectively. The full listing is available 
electronically. 

4. An old and massive bulge? 

Before discussing the stellar populations (Sect. 14. Il l and dynam- 
ics (Sect. 14. 2t of the central 300 arcsec of M31, it is important 
to assess what is the light contribution of the bulge and the disk 
of M31 as a function of distance and position angle. According 
to the decomposition of Kormendy and Bender (1999 ), based on 
a Sersic plus exponential law fit to the MJ photometry and very 
similar to the one of Kent (I1989bl l. the disk makes 1 % of the to- 



tal light at the center, 10% at 100 arcsec and 31% at a distance of 
300 arcsec along the major axis, and 9% and 25% along the mi- 
nor axis respectively. The disk light amounts to j 55% of the total 
at the end of the MJEE slit, 500 arcsec from the center. The re- 
sults discussed in Sect. l4.1l show that generally old stellar popu- 
lations are found out to 300 arcsec from the center. In agreement 
with Stephens et al. (120031 1. who start detecting (younger) disk 
stars only at distances larger than 7 arcmin from the center, we 
see the effect of the increasing importance of the younger disk 
component only at distances of 400-500 arcsec from the center. 
In contrast, the kinematic imprint of the disk stars is probably 
visible already at distances > 100 arcsec, especially on the ma- 
jor axis (see Fig.|3]and Sect. |3]l. We discuss it in Sect. 14.21 
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Fig. 3. Continued 

Table 2. Format example of the measured stellar kinematics as a function of distance from the center (positive: east, negative: west) 
for the different position angles. The full table is available electronically. 



R PA KiMrs stars ^^4 

_n (km/s) (km/s) 

-289 MJ -404.14 ±2.21 143.23 ± 2.60 -0.061 + 0.014 -0.002 ± 0.014 



Table 3. Format example of the measured gas kinematics as a function of distance from the center (positive: east, negative: west) 
for the different position angles. The full table is available electronically. 



R PA Vg„s cTgcs up [Om]/UJ3 [NI]/H/3 

(") (km/s) (km/s) A 

■^289 MJ -458.7 + 4.9 94 ±5.6 0.33 ± 0.046 2.61 +0.15 0.01 ± 15.7 



4. 1 . Stellar populations 

We study the stellar populations of the inner 5 arcmin of M3 1 
(and out to 8 arcmin along the major axis) using the simple 
stellar population (SSP) models of Maraston ([19981 12005 ) with 
a Kroupa (I2001I I IMF, and the modeling of the Lick indices 
(LSSP) with a-elements overabundance of Thomas, Maraston 
and Bender (2003). The models cover ages t up 15 Gyr, metallic- 
ities [Z/H] from -2.25 dex to -1-0.67, and overabundances [a/Fe] 
from -0.3 to H-0.5. We spline-interpolate the LSSP models on a 



finer grid in f, [Z/H] and [a/Fe] and determine values of f, [Z/H] 
and [c/Fe] that minimize the;^'(r)^: 



X(rf = AHpirf + AMgbirf + AFe5015(rf+ 

AFe521Q(rf + AFe5335{rf + AFe5'\Q6{rf+ 
AFe51Q9{rf + AFe5782(rf 

at each radius r, where Alndexir)^ is: 



AIndex{r) 



(Index(r) — mlndex(t, [Z/H], [a/Fe])\ 
dindexir) / 



(1) 



(2) 



R.P. Saglia et al.: The old and heavy bulge of M31 



t/2 



a 



b 



> 
xn 



m 



b 



m 

CO 




10 100 

r [arcsec] 



10 100 
r [arcsec] 



10 100 
r [arcsec] 



Fig. 5. The gas kinematics along different position angles, folded with respect to the center, antisymmetrically for recessional 
velocities, symmetrically for velocity dispersions. Only data with errors smaller than 50 km/s are shown. Filled circles refer to 
radii with increasing RA, stars to radii with decreasing RA, with the same convention and zero-point shifts on radial velocities of 
Fig. [3] The vertical dotted lines mark the transition from the central slit data to those from the outer two. For the major axis, a 
further vertical dotted line marks the outer MJEE dataset. The (red) short-dashed vertical lines mark the seeing values of the central 
slits (see Table [T]). The data of Pellet IT976l l are also shown for comparison as open circles (for increasing RA) and triangles (for 
decreasing RA). His 28 A data are plotted green, the 135 A red. We show his PA=45 data as MJ, his PA=128 as MN, his PA=68 as 
P30. Moreover, the data at PA=38 (his major axis) are shown blue also as MJ. 

Table 4. Format example of the measured Lick indices as a function of distance from the center (positive: east, negative: west) for 
the different position angles. The full table is available electronically. 



R PA 


W 


Mgb 


Fe5015 


Fe5270 


Fe5335 


Fe5406 


Fe5709 


Fe5782 


(") 


A 


A 


A 


A 


A 


A 


A 


A 


-289 MJ 


1.65 ± 0.21 


4.04 ±0.17 


5.07 ± 0.40 


2.94 ±0.14 


2.50 ±0.32 


1.53 ±0.12 


0.78 ±0.11 


0.78 ± 0.08 



where dlndex{r) is the error on the index lndex{r) measured at 
distance r and mlndex(t, [Z/H], [a/Fe]) is the value predicted 
by the LSSP models for the given set of age, metallicity and 
overabundance. Note that we do not extrapolate to values out- 
side the model grid. Therefore the minimum solutions are 
sometimes found at the edges of the parameter space. Once the 
set (tmi„(r), [Z/H],„i„(r), [a/Fe],„i„(r)) that minimizes x(rf has 



been determined, we compute the values of the colors and mass- 
to-light (M/L) ratios of the corresponding SSP models. Errors on 
all quantities are derived by considering the minimal and maxi- 
mal parameter variations compatible with Ax^ - - Xmi„ - 1 ■ 

Fig. [TT] shows the resulting age, metallicity, overabundance, 
color and M/L profiles along the different slit positions. The 
corresponding best-fit lines to the measured indices profiles are 
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Fig. 7. The comparison between the line strengths measured here and the Lick values for 20 stars. The dotted line shows the one-to- 
one correlation, the full line the best fitting relation ULick) - ax Lick(Us) + b. 



shown in Fig.|9] We fit 8 Lick indices to derive three parameters, 
so the expected should be around 5. This matches rather well 
the derived values (see bottom plot of Fig.|9l), indicating that the 
SSP models are a good representation of the data within their 
sizeable errors. 

The six age, metallicity and overabundance slit profiles are 
very similar when overplotted as a function of radius. This is a 
consequence of the approximately constancy of the indices along 



circles discussed above and indicates that the region with homo- 
geneous stellar populations is not as flattened as a disk compo- 
nent would be, in agreement with what found by Davidge ( 19971 ) 
for the inner 30 arcsec. Fig. [12] shows the age, metallicity, over- 
abundance, and MILr profiles averaged over position angles and 
binned logarithmically in radius (for the MN profile, only the 
central slit data are included in the average). 
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Fig. 9. The measured Lick indices profiles along the 6 slit positions of M3 1 . The vertical dotted lines mark the ends of the central 
slits. For the major axis, a further vertical dotted line marks the outer MJEE dataset. The lines show the best-fit LSSP models (see 
discussion in Sect. 14.1b . The bottom plots show the respective \ogx^, coded by color, and the total \ogx^ in black. The long-dashed 
line shows the expected value for the given degrees of freedom (5). The open symbols in P30 show the data of Davidge (119971) . 



The picture emerging from Fig. [TT| and [12] is the following. 
At radii larger than a few arcsec and out to a; 5 arcmin from the 
center, the stars of the bulge of M3 1 are on average almost as 
old as the universe, tend to solar metallicity and are a overabun- 
dant by 0.2-0.25 dex. This is in agreement with Stephens et al. 
(I2003J . who find that the near infrared luminosity function of the 



red giant stars of the bulge of M3 1 is indistinguishable from the 
one of the Milky Way bulge. Of course, the point to point scat- 
ter in age in the single slit profiles is large and formally the 1- 
sigma lines would allow sometimes ages as small as 5 Gyr. This 
is driven by the scatter in the H/?, that ultimately correlates with 
the emission correction applied. There is a metalhcity gradient 
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Fig. 9. Continued. 



of 0.2 dex per radial decade out to 40 arcsec, and essentially no 
gradient at larger distances, where solar metallicity is reached. 
There is no obvious radial variation in age or overabundance. 
This translates in moderate u-g and M/Lb gradients, and almost 
no detectable gradients in the g-r, r-i, M/Ly and M/Lr. The exact 
values of the stellar mass-to-light ratios depend on the assumed 
IMF (Kroupa 2001 ), that has a slope a = 1.3 for the mass range 
0.1 < M < O.SMq and a = 2.3 at larger masses. Changing the 
break mass from Q.5Mq to IMq lowers the M/Lr from 4 to 3. 



Along the major axis and at distances larger than ^ 6 arcmin, 
the mean age of the stellar population drops to < 8 Gyr, again in 
agreement with Stephens et al. (2003), who find disk stars only 
at distances larger than 7 arcmin from the center The metallicity 
is slightly supersolar (by ^ -i-O.l dex). The [a/Fe] overabunce 
remains high. As a consequence, the derived mass-to-light ratio 
drops to M/Lr < 3Mq/Lq. 

The situation is different in the inner arcsecs, where the stel- 
lar population age drops down to 4-8 Gyr, the metallicity in- 
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Fig. 10. The EW of the H/3 emission line in A and the ratios [OIII]/Hj6 and [NI]/Hj0 of the emission EWs in logarithmic units as a 
function of radius along different position angles. The large filled circles show the points with signal to noise ratio larger than 4. 
The vertical dotted lines mark the ends of the central slits. For the major axis, a further vertical dotted line marks the outer MJEE 
dataset. 



creases to above 3 times solar, but the overabundance remains 
at » 0.2 dex. This is in broad agreement with the conclusions of 
Davidge ( 1997), who however hints to a possible further increase 
of the overabundance towards the center, not seen here. 



We compare the predicted Sloan color profiles to the ones de- 
rived from the Sloan survey, that have been dereddened using the 
extinction map discussed in Montalto et al. ("2009). Following 
Maraston (private communication) we apply a correction of - 
0.05 mag to the model g - r and of +0.07 mag to the model r-i 
colors. The comparison is surprisingly good in u-g and good in 
r-i, with deviations < 0.1 mag. Some residual systematic devia- 
tions are still present in g-r (where the models are too red) and r-i 
(where the models are slightly too blue), in qualitative agreement 
with the discussion of Maraston et al. (2009). This overall match 
provides an independent assessment of the line index analysis. 



4.2. Axisymmetric dynamical modeling 

The kinematics shown in Fig. |3]describe the light averaged stel- 
lar motions of the bulge and disk components. In the inner ar- 
cmin this follows closely the bulge kinematics, but at larger dis- 
tances the contamination by the disk becomes more and more 
important. Therefore, a dynamical model for the sum of the two 
components should be compared to Fig. [3] Qualitatively, we can 
expect the stellar disk to rotate much faster, possibly up to the 
150 km/s observed in the gas (see Fig.|5]l, and to be colder (with 
velocity dispersions of the order of 60-70 km/s) than the bulge. 
As a consequence, the true rotational velocity of the bulge could 
be smaller than, and the bulge velocity dispersion could be larger 
than shown in Fig.|3] Moreover, the line of sight velocity distri- 
bution in the outer regions will be dominated in the high velocity 
tail by disk stars. This can balance the intrinsic asymmetry of the 
bulge (that anticorrelates with the mean rotational velocity), pos- 
sibly explaining the change in sign of the Hj, coefficient observed 
near the major axis at distances larger than 100 arcsec. We have 
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Fig. 11. The age, metallicity, overabundance, Sloan colors u-g (blue), g-r (green), r-i (red) and Johnson M/L profiles in the B (blue), 
V (green) and R (red) bands along the different slit positions. The dotted lines show the 1-sigma errors. The vertical dotted lines 
mark the ends of the central slits. For the major axis, a further vertical dotted line marks the outer MJEE dataset. The measured 
color profiles, corrected by extinction (Montalto et al. l2009l l are shown as crosses. 



tested that this scheme could work with a simple bulge plus disk 
kinematical model that follows the light decomposition of Kent 
( |1989bl l, but we postpone a quantitative analysis to Morganti et 
al. (l2009] l. 

The blue lines in Fig. [3] show the fits to the data of McElroy 
( fT983D by Widrow et al. (|2003l l, model A (full) and model Kl 



(dashed). From the discussion given above, it is clear that it 
is dangerous to assume that McElroy's data describe the kine- 
matics of the bulge component alone, as done by Widrow et al. 
(I2003I I. Indeed, McElroy (|19831 l recognises that he might be un- 
derestimating the true stellar velocity dispersion due to how his 
analysis method reacts to disk light contamination (see his Table 
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Fig. 11. Continued. 



3), but concludes that the effect is not large on the basis of the 
disk to bulge decomposition available at the time, that severely 
downplayed the role of the disk in the inner 5 arcmin. 

Letting aside this intrinsic problem, the best-fit model A of 
Widrow et al. (I2003I I predicts a mass-to-light ratio for the bulge 
(M/Lr = 2.7Mq/Lq) that is smaller than the one of the disk 
(M/Lr = 4AMq/Lq), and this despite overestimating the ve- 
locity dispersions along the minor axis. Only model Kl, well 
above McElroy crs, delivers M/Lr = 4.SMq/Lq for the bulge. 



Widrow et al. go on discussing why it is so, not taking into ac- 
count the possibility that McErloy velocities are in fact systemat- 
ically too small. Of course Widrow et al. model A fails to match 
our kinematics, while model Kl slightly overestimates our ve- 
locity dispersions. Fig.fTTIshows that for r > 10 arcsec the anal- 
ysis of the Lick indices predicts M/Lr ~ 4 - 4.5Mq/Lq, in line 
with model Kl. As a consequence, the microlensing events es- 
timated by Widrow et al. (2003), that are based on their model 
A, are clearly underestimates. Similarly, an upward revision of 
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Fig. 4. Top: the full line shows the velocity dispersion result- 
ing by averaging all measured cr within the isophote with semi- 
major axis a in arcsec, with weights equal to the surface luminos- 
ity at each point. The dotted line shows the same for the quan- 
tity Vcr^ + V^, where V is the measured line of sight velocity. 
Bottom: the same computed by weighting according to the total 
isophotal light. 



the event rate estimates of Riffeser et al. (12006') is needed, since 
this paper assumes a stellar mass-to-light ratio for the bulge 
of M/Lr = 2.6. The extremely low M/Lr value (0.8 M^/Lq) 
favoured by Chemin et al. (2009 ) based on a simple mass model- 
ing of the HI velocities is completely ruled out and points to the 
inadequacy of underlying modeling assumptions (circular mo- 
tions and spherical mass distributions). 

Finally, a word on the triaxiality of the bulge or the possi- 
ble presence of a bar. A rotating axisymmetric structure shows 
equal velocity curves along slits inclined at positive or negative 
angles with respect to the major or minor axis. In general, this 
does not happen if the structure is triaxial or a bar is present 
(see Fig. 12 of Athanassoula & Beaton |2006l l. In the case of the 
bulge of M3 1 the maximal rotation velocity achieved on the slits 
M60 and M30 are larger than the ones reached on slits P60 and 
P30 respectively, or, in other words, the kinematic minor axis of 
the M31 bulge does not coincide with its isophotal minor axis. 
We quantify this effect by measuring how the position angle of 
the kinematic major axis varies as a function of semi-major dis- 
tance. At each semi-major distance a from the center we find the 
parameters Vmaxio) and (poia) that minimize thexi'^)^ function: 
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y,(fl) - V„,ax{a) cos(0i - <f>o{a)) 



6Vi{a) 



(3) 



where the twelve angles 0, correspond to the ones of the slits 
shown in Fig. [1] and listed in Table [1] and y,(fl) are the stel- 
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Fig. 6. The variation of the molecular Mgl and Mg2 (in mag) 
and atomic (in A, all the rest) line strengths as a function of 
position on the slit measured for the star HD72324 (= HR 3369). 
The dotted line shows the values of Worthey et al. (il994j . 



lar velocities relative to the center shown in Fig. [3] interpolated 
at the distance proper to the isophote ellipse determined below 
(see top panel of Fig. [T3] ). In Eq. [3] we exclude the velocities 
measured along the minor axis at distances r > 100 arcsec, 
where the slits are shifted with respect from the center A per- 
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Fig. 8. The effect of the gas emission correction on the H/3, 
Fe5015 and Mgb indices along the major axis of M31. The plot 
shows the differences A/ = Icor - haw between the emission cor- 
rected index Icor and the values haw measured before the correc- 
tion. The full lines show the EW of the measured emissions in 
Hj6 and [NI] and the difference between the EWs of the redder 
and bluer emission lines of the [OIII] doublet (see text). 



fectly axisymmetric system has 0o = PAmj, where PAmj is the 
position angle of the isophotal major axis. In general, the func- 
tion Vmax cos(0,- - 0()) is a fair representation of the measured 
data, with rms deviations of 8 km/s. Fig. [13] bottom, shows the 
resulting profile <poia), together with PAMj{a) as derived from 
a 2MASS K image of M31 using the ellipse fitting package 
of Bender and Mollenhoff ( |1987i ). The corresponding ellipticity 
profile is shown at the top. The isophotal PA and ellipticity pro- 
files agree with Beaton et al. ( I2007I I. Both the isophotal and the 
kinematic PAs move with distance from the center. But, while 
the isophotal PA oscillates around the value of 48° (the major 
axis PA of the M31 bulge, see Table [1]), the kinematic PA drifts 
away from it, approaching 32°, 6 degrees smaller than the PA of 
the major axis of the disk. This suggests that the observed drift is 
due only partially to the increasing importance of the stellar disk 
rotation and might point to an intrinsic triaxiality of the bulge 
and/or the presence of a bar. 

However, only a proper dynamical model that takes into ac- 
count the disk contribution, on top of the triaxiality of M3 1 's 
bulge and its possible figure rotation, and fits not only the light 
distribution, mean velocities and velocity dispersion fields, but 
also the higher moments of the LOSVDs, can address these is- 
sues quantitatively. We postpone this discussion to Morganti et 
al. (l2U09b . 
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Fig. 12. The age, metallicity, overabundance, and M/Lr profiles 
(full lines) averaged over position angles and binned logarith- 
mically in radius. The dotted lines and the bars show the corre- 
sponding rms and errors on the mean, respectively. 



4.3. The sources of ionizing radiation 



Sarzi et al. ( I2009I I discuss in detail how to exploit the [OIII]/HyS 
vs. [NI]/HjS diagram to constrain the properties of the ionizing 
sources for the gas in elliptical and lenticular galaxies. Their 
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Fig. 13. Top: the ellipticity profile as a function of the semi- 
major axis distance a. Bottom: the position angle of the kine- 
matic major axis (poia) (crosses) and of the isophotes (filled 
circles) as a function of a. Both the ellipticity profile and the 
isophotes PAs are derived from a K band 2MASS image of M3 1 . 



Fig. 1 shows that the diagram neatly separates the regions where 
dusty AGNs, shocks and startburst activity are responsible for 
the production of emission lines, in analogy with the stan- 
dard [OI]/Ha vs. [011I]/Hj6 diagram of VeiUeux and Osterbrook 

(fT9871 l. 

As discussed by Sarzi et al. (2009V the constancy of the Wfi 
EW with radius indicates that the Hfi emission flux follows the 
stellar light distribution. To first approximation, this should ap- 
ply also to the case discussed here. Values of log [0111]/Hy6 * 0.5 
with low [Nl]/Hy6 ratios occupy the region where shocks are re- 
sponsible for the gas excitation (see Fig. 1 of Sarzi et al. 120091 ). 
In contrast, [OIII]/HjS ratios as high as 10, like the ones measured 
in the inner arcsecs of M31, fall in the region where photoion- 
ization by a central AGN is the working mechanism. 



5. Discussion and Conclusions 

We have presented new optical spectroscopic observations of the 
bulge of M3 1 . We have measured the stellar and gas kinematics, 
emission line strength ratios and absorption Lick indices profiles 
along 6 position angles out to distances of 5 arcmin from the 
center Along the major axis we probed regions out to 8 arcmin. 
We show that the old kinematics of McElroy ( 119831 ) provides 
too small velocity dispersions (by up to 30%), therefore biasing 
the dynamical modeling towards too small bulge stellar masses 
(by a factor 2). Moreover, the new higher averaged velocity dis- 
persion predicts a mass for the central supermassive black hole 
of M3 1 that is only a factor 2 below what measured. The new 
velocity dispersion profiles are now in better agreement with ax- 
isymmetric dynamical models with large bulge mass-to-light ra- 
tio (Widrow et al. |2003l l, that now match the values derived from 
stellar population models {M/Lr ~ 4 - 4.5Mq/Lq, see below). 
This implies an upward revision of the predicted self-lensing mi- 
crolensing event rate of Widrow et al. (2003 ), and Riffeser et al. 
( 120061) . that are based on lower stellar mass-to-light ratios. 

The inner (r < 100 arcsec) bulge is slowly rotating, with a 
V/cr 0.2. At distances from the center larger than ^ 100 arcsec 
the measured kinematics becomes increasingly influenced by the 
rapidly rotating stellar disk. Therefore, the observed variation of 
the kinematic position angle is suggestive of bulge triaxiality, but 
needs a proper dynamical modeling of both disk and bulge com- 
ponents to be quantified. The measured gas kinematics confirms 
the well studied large scale disk rotation. However, a more com- 
plex structure, with gas minor axis counter-rotation, is detected 
in the inner bulge. This might be evidence for a (recent) minor 
merger, possibly connected to the younger stellar population de- 
tected in the inner arcsecs of the galaxy discussed below. 

The analysis of eight Lick index profiles shows that the bulge 
of M3 1 is old, of solar-metallicity and a factor 2 overabundant 
in a-elements, in agreement with studies of its resolved stellar 
component (Stephens et al. 2003). The line indices and stel- 
lar population parameters appear approximately constant on cir- 
cles, i.e. their isocontours are rounder than the galaxy isophotes, 
as seen in many ellipticals and bulges (Kuntschner et al. 120061 
Falcon-Barroso et al. 120061) . Together with the derived old ages, 
this confirms that the stellar disk out to 5 arcmin from the cen- 
ter is either old (similar to what found for other spiral galax- 
ies, Peletier and Balcells 1996 ) or not sufficiently probed by the 
spectral features considered here. However, we do detect smaller 
ages (< 8 Gyr) along the major axis at distances > 6 arcmin. The 
u-g SLOAN colors predicted from our stellar population analysis 
match reasonably well with the observed ones. The redder col- 
ors g-r and r-i are systematically offset, a well known problem of 
the flux calibration of current simple stellar population models 
(Maraston et al. |2009l l. The derived mass-to-light ratios (in the 
Johnson R band and with a Kroupa IMF we get M/L =s 4 - 4.5) 
agree with the dynamical estimates (see above). They drop to 
< 3Mq/Lq along the major axis at distances > 6 arcmin, where 
the disk light starts to dominate. 

In the inner arcsecs the situation changes and a population 
with a light-weighted younger age 8 Gyr inside the see- 
ing disk or 2 arcsec, with values as low as 4 Gyr) and metal 
richer (^ 3 times solar) appears. This agrees with the findings 
of Davidge ( 119971 1. In addition, the emission line EW ratios 
[OIII]/Hj0 increase in this region. Their values are compatible 
with being excited by shocks in the main body of the bulge, but 
near the center they increase to levels pointing to the presence 
of an AGN-like photo-ionizing source. Combined with the de- 
tection of counterrotating gas along the minor axis of the galaxy 
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(see above), this suggests that a gas-rich minor merger probably 
happened some 100 Myr ago, that triggered an episode of star 
formation and possibly boosted the nuclear activity of the central 
supermassive black hole of M3 1 . We estimate how broad a range 
of star burst ages and masses can be to give the measured mean 
value of 8 Gyr in the inner 2 arcsec, when superimposed to the 
old bulge stars background. To this purpose we compute com- 
posite spectra of an old (12.6 Gyr) plus a young (from 100 Myr 
to 4 Gyr) simple stellar population, using the Vazdekis d 19991 1 
library, and measure their SSP ages through the analysis of the 
Lick indices observed here. We find that global ages smaller than 
8 Gyr are found when considering a young component younger 
than ^ 600 Myr and a mass fraction lower than 10%. Higher 
mass fractions are possible for older ages. From Kormendy and 
Bender ( 11999 ), inside an aperture of 2 arcsec diameter we mea- 
sure a V mag of 12.5, or 5 x lO^L©. In this region we estimate 
MjLv ~ 4Mq/Lq and therefore an enclosed mass of 2 x lO^M©. 
As a consequence, x lO^M© of some 100 Myr old stars would 
be needed. Note that in the inner nucleus of M31, at fractions 
of an arcsec, a disk of 200 Myr old stars is found (Bender et al. 
I2005I I. with a mass of « 4200Mo. Of course the spatial resolu- 
tion of our spectra is too low to probe this scale. Moreover, our 
result, based on the H/3 line as an age tracer, depends heavily 
on the details of the emission correction and might be affected 
by the bad column's interpolation (see Sect. |2]i. Spectra of the 
higher order Balmer lines, taken with better seeing, are needed 
to improve our conclusions. 

In a second paper (Morganti et al. 120091), a dynamical model 
of the data presented here, that takes into account the contribu- 
tions of the bulge and disk components, will assess in a quanti- 
tative way the bulge triaxiality issue and will give a proper esti- 
mate of the microlensing event rates due to self-lensing. 
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